##### SCRIPT QUE GERA A FIGURA 4 #####

# 1. CARREGA PACOTES ------------------------------------------------------

library(tidyverse)
library(janitor)
library(naniar)
library(mice)
library(xlsx)

# 2. CARREGA DADOS ------------------------------------------------------

#' Nesta etapa, foi implementada a leitura dos dados combinados de Enade-RAIS (filtrado obs com salario)

# 3. IMPLEMENTA MODELOS DE REGRESSAO OLS --------------------------------

# 3.1 Modelo separado por área ----

imp_med <- filter(dados_graduados_empregados, area_curso_agreg == "Medicina")
imp_soc <- filter(dados_graduados_empregados, area_curso_agreg == "Ciências Sociais Aplicadas")
imp_ctem <- filter(dados_graduados_empregados, area_curso_agreg == "CTEM")

imp_dir <- filter(dados_graduados_empregados, area_curso_agreg == "Direito")
imp_edu <- filter(dados_graduados_empregados, area_curso_agreg == "Educação")
imp_eng <- filter(dados_graduados_empregados, area_curso_agreg == "Engenharia")

imp_hum <- filter(dados_graduados_empregados, area_curso_agreg == "Humanidades")
imp_sau <- filter(dados_graduados_empregados, area_curso_agreg == "Saúde e Bem-Estar")
imp_tec <- filter(dados_graduados_empregados, area_curso_agreg == "Tecnólogos")

f <- "log(salario_hora) ~ tp_sexo + raca_cor + idade_faixa + edu_fam + nt_ce_quartil_area + nt_fg_quartil + trab_grad + igc_faixa_4c + setor_ies_bin + regiao_curso"

model_med <- with(imp_med, lm(as.formula(f)))
model_soc <- with(imp_soc, lm(as.formula(f)))
model_ctem <- with(imp_ctem, lm(as.formula(f)))
model_dir <- with(imp_dir, lm(as.formula(f)))
model_edu <- with(imp_edu, lm(as.formula(f)))
model_eng <- with(imp_eng, lm(as.formula(f)))
model_hum <- with(imp_hum, lm(as.formula(f)))
model_sau <- with(imp_sau, lm(as.formula(f)))
model_tec <- with(imp_tec, lm(as.formula(f)))

combine_med <- pool(model_med)
combine_soc <- pool(model_soc)
combine_ctem <- pool(model_ctem)
combine_dir <- pool(model_dir)
combine_edu <- pool(model_edu)
combine_eng <- pool(model_eng)
combine_hum <- pool(model_hum)
combine_sau <- pool(model_sau)
combine_tec <- pool(model_tec)

# 4. CALCULA SALARIOS PREDITOS -------------------------------------

# 4.1 Por CE, area e setor da IES -----

# Funcao que cria dados-base para calculo de probabilidades preditas

prepara_new_data <- function(x = data, y) {
      x %>% 
            dplyr::distinct(area_curso_agreg) %>%
            filter(area_curso_agreg == y) %>% 
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            dplyr::arrange(area_curso_agreg) %>% 
            dplyr::mutate(tp_sexo = "Feminino",
                          raca_cor = "Negro/Indígena",
                          edu_fam = "Ensino Médio",
                          nt_ce_quartil_area = rep(c("q1", "q4"), 2),
                          nt_fg_quartil = "q3",
                          trab_grad = c("Não trabalha"),
                          igc_faixa_4c = "3",
                          idade_faixa = "25-",
                          setor_ies_bin = rep(c("Pública", "Privada"), each = 2),
                          regiao_curso = "SE") %>% 
            mutate_if(is.character, as.factor)
}

new_data_soc <- prepara_new_data(dados_graduados_empregados, "Ciências Sociais Aplicadas")
new_data_cte <- prepara_new_data(dados_graduados_empregados, "CTEM")
new_data_dir <- prepara_new_data(dados_graduados_empregados, "Direito")
new_data_edu <- prepara_new_data(dados_graduados_empregados, "Educação")
new_data_eng <- prepara_new_data(dados_graduados_empregados, "Engenharia")
new_data_hum <- prepara_new_data(dados_graduados_empregados, "Humanidades")
new_data_med <- prepara_new_data(dados_graduados_empregados, "Medicina")
new_data_sau <- prepara_new_data(dados_graduados_empregados, "Saúde e Bem-Estar")
new_data_tec <- prepara_new_data(dados_graduados_empregados, "Tecnólogos")

# Funcao que calcula probabilidades preditas

predict_prob <- function(x, y, z) {
      p <- x$analyses[[1]]
      p$coefficients = summary(y)$estimate
      predict(p, newdata = z,
              type = "response", se.fit = TRUE)
}

pred_soc <- predict_prob(model_soc, combine_soc, new_data_soc)
pred_cte <- predict_prob(model_ctem, combine_ctem, new_data_cte)
pred_dir <- predict_prob(model_dir, combine_dir, new_data_dir)
pred_edu <- predict_prob(model_edu, combine_edu, new_data_edu)
pred_eng <- predict_prob(model_eng, combine_eng, new_data_eng)
pred_hum <- predict_prob(model_hum, combine_hum, new_data_hum)
pred_med <- predict_prob(model_med, combine_med, new_data_med)
pred_sau <- predict_prob(model_sau, combine_sau, new_data_sau)
pred_tec <- predict_prob(model_tec, combine_tec, new_data_tec)

# Funcao que gera tabela com probabilidades preditas

prepara_tabela_2 <- function(x, y, z) {
      data.frame(z[,c("area_curso_agreg", "nt_ce_quartil_area", "setor_ies_bin")],
                 pred = x) %>%
            mutate(pred.fit = exp(pred.fit),
                   pred.se.fit = exp(pred.se.fit)) %>% 
            mutate(area_curso_agreg = y) %>% 
            select(1, 2, 3, 4, 5)
}

tab_ce <- bind_rows(
      prepara_tabela_2(pred_soc, "SOC", new_data_soc),
      prepara_tabela_2(pred_cte, "CTEM", new_data_cte),
      prepara_tabela_2(pred_dir, "DIR", new_data_dir),
      prepara_tabela_2(pred_edu, "EDU", new_data_edu),
      prepara_tabela_2(pred_eng, "ENG", new_data_eng),
      prepara_tabela_2(pred_hum, "HUM", new_data_hum),
      prepara_tabela_2(pred_med, "MED", new_data_med),
      prepara_tabela_2(pred_sau, "SAU", new_data_sau),
      prepara_tabela_2(pred_tec, "TEC", new_data_tec))

# 5. GERA GRAFICO --------------------------------------------------------

# Cria grafico

fig_4 <- tab_ce %>% 
      unite(setor_quartil, setor_ies_bin, nt_ce_quartil_area) %>% 
      ggplot(aes(area_curso_agreg, pred.fit, fill = setor_quartil)) +
      geom_bar(stat = "identity", position = position_dodge(.7)) +
      scale_y_continuous(breaks = seq(0, 40, by = 5)) +
      geom_errorbar(aes(ymin=pred.fit-pred.se.fit, ymax=pred.fit+pred.se.fit), width=.2,
                    position=position_dodge(.7)) +
      theme_bw() +
      scale_fill_grey(start = 0.8, end = 0.3) +
      labs(x = "", y = "Salário-hora predito") +
      theme(legend.position = "bottom", 
            legend.title = element_blank(),
            text = element_text(size = 18, color = "black"),
            axis.text = element_text(size = 18, color = "black"),
            axis.title = element_text(size = 18, color = "black"))

# 6. EXPORTA GRAFICO -----------------------------------------------------

jpeg(filename = "output/Figura 4.jpeg", width = 582, height = 388)
fig_5
dev.off()
